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Abstract 



An important issue with oversampled FIR analysis filter banks (FBs) is to determine inverse synthesis 
■ FBs, when they exist. Given any complex oversampled FIR analysis FB, we first provide an algorithm 

to determine whether there exists an inverse FIR synthesis system. We also provide a method to ensure 
the Hermitian symmetry property on the synthesis side, which is serviceable to processing real-valued 



signals. As an invertible analysis scheme corresponds to a redundant decomposition, there is no unique 
inverse FB. Given a particular solution, we parameterize the whole family of inverses through a null 
space projection. The resulting reduced parameter set simplifies design procedures, since the perfect 
reconstruction constrained optimization problem is recast as an unconstrained optimization problem. The 

in 

\Q . design of optimized synthesis FBs based on time or frequency localization criteria is then investigated, 

m 



using a simple yet efficient gradient algorithm. 



Index Terms 



Oversampled filter banks, inversion, filter design, optimization, time localization, frequency localiza- 

, 

■ tion, lapped transforms, modulated filter banks. 



I. Introduction 

Since the 70s, filter banks (FBs) have become a central tool in signal/image processing and commu- 
nications: lapped or discrete wavelet transforms can be viewed as instances of FB structures. Likewise, 
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oversampled FBs (OFBs) constitute an extensively studied instance with remaining open questions. Their 
development came along under a variety of different appellations, to name a few: general Analysis- 
Synthesis Systems 12, DFT (discrete Fourier transform) with stack-shift capability, Overlap-Add or 
Generalized DFT, underdecimated systems, oversampled harmonic modulated filter banks flU, Q, com- 
plex lapped transforms (4J, generalized lapped pseudo-biorthogonal transform, etc. 

In a more generic form, OFBs have received a considerable attention both theoretically and in many 
applications, in the past ten years, following their association with specific types of frames (H, ||5], 0. 
Their design flexibility, improved frequency selectivity and increased robustness to noise and aliasing 
distortions have made them useful for subband adaptive filtering in audio processing (7), noise shaping 
flU, denoising Q, multiple description coding (9[, echo cancellation flOl . multiple antenna code design 
ifTll . channel equalization El, (BJ, HH or channel coding |fT5Tl . 

Two major problems arise when resorting to OFBs: (J) the existence of an inverse for the analysis OFB 
achieving perfect reconstruction (PR) and («) the determination of an "optimal" synthesis FB. Since the 
additional degrees of freedom gained through redundancy may increase the design complexity, several 
works have focused on FBs modulated with a single lfT6*l . ifTTl or multiple windows lfl8l . More general 
formulations are based on factorizations of OFB polyphase representations with additional constraints 
(restricted oversampling ratios, symmetry, realness or filter length) into a lattice |[T9l , EOl , BT1 . E2l or a 
lifting structure ll23l . Constructions with near perfect reconstruction (relaxing the PR property) have also 
been proposed EH, [25]], HQl, EH. In (271, EH, EH, more involved algebraic tools (such as Grobner 
bases) have also been employed. Recently, Chai et al. have proposed a design based on FB state-space 
representations ll30l . The design may use different kinds of optimization criteria based on filter regularity 
or more traditional cost functions based on filter shape (subband attenuation ED, iTTOl . coding gain ||3T1 ). 
Most of those synthesis FB designs rely on minimum-norm solutions. An interesting approach combining 
the latters with a null space method was successfully pursued by Mansour ll32l for synthesis window 
shape optimization in a modulated DFT FB. 

Within the compass of the proposed work is a relatively generic construction and optimization of 
oversampled synthesis filter banks with Finite Impulse Response (FIR) properties at both the analysis 
and synthesis sides. We can additionally impose a practically useful Hermitian symmetry on the synthesis 
side. This work extends the results given in two previous conference papers |[33l , ||34l . A special case 
has been judiciously devised in 11221 . for specific filter length and redundancy factor allowing closed 
form expressions for two design criteria. In Section [TT] we recall the polyphase notation used throughout 
this paper. Given arbitrary FIR complex oversampled analysis FB, we first describe in Section MI- A I 
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a simple algorithm to test whether it is FIR invertible or not, based on known results on polynomial 
matrices ll35l . ll36l . The standard Moore-Penrose pseudo-inverse (PI) solution ll37l is studied in Section 
IIII-B I In Section IIII-CI a method is supplied to enforce an Hermitian symmetric FB, which is useful 
for real data analysis, processing and synthesis. In Section |IVl the problem of the optimal design of the 
synthesis FB is addressed. Although optimization can be studied both on the analysis and synthesis sides 
ll38l . |39l , we consider here a given analysis FB and work on the synthesis side. We derive in Section 
IIV-AI an efficient parameter set size reduction for this purpose. Using time or frequency localization 
criteria, we then reformulate in Section IIV-BI the constrained optimization problem as an unconstrained 
one for both the general and Hermitian symmetric cases. After describing the optimization process, 
we illustrate, in Section |V] the different methods proposed for the inversion and optimization on three 
classical oversampled real and complex FB types. 



A. Notations 



II. Problem statement 
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Fig. 1. Oversampled M -channel filter bank. 



Lapped transforms ||40l were introduced in HTI to avoid blocking artifacts in audio processing. 
Similarly for images, they reduce tiling effects produced by classical block transforms (as can be seen 
in the JPEG image compression format). Lapped transforms belong to the class of FBs, such as the one 
represented in Figure [1] with a decimation factor N smaller than the length of each filter. The filters, 
whose impulse responses are denoted by (/ij)o<i<M> are supposed of finite length kN with k an integer 
greater than or equal to 2. We therefore consider k overlapping blocks of size N. 

A signal (x(n)) ne z is decomposed by M filters; since the decimation factor is N, the overall redun- 
dancy of the transform is M/N = k'. In this paper, we investigate the oversampled case, i.e. k' > 1. 
The M outputs of the analysis FB are denoted by (v7i(n))o<i<ju- With these notations, the outputs of 
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the analysis FB are expressed, for all i G {0, . . . , M — 1} and n G Z, as 

N-l 

y i {n) = Y J h t {p)x{Nn-p) = Y J Y. hi(N£ + j)x(N(n - 1) - j). (1) 

V i j=0 

B. Polyphase formulation 
LetH(£) = [hi(N£ + j)} 

o<i<M,o<j<N> ^ G {0, . . . , A; — 1} be the k polyphase matrices obtained from 
the impulse responses of the analysis filters. We also define the polyphase vector from the input signal 
x{n): Vn E Z, x(n) = (x(Nn — j))o<j<N> l eacrm g to concisely rewriting (fl]) into a convolutive form: 

y(n) = (y (n), . . . , y M -i(n)) T = £ H(*)x(n = (n), (2) 

£=0 

where T is the transpose operator. Thus, © can be reexpressed as: y[z] = H[z]x[z], where H[z] = 
Y^.ZoH(l)z~^ is the M x N polyphase transfer matrix of the analysis FB and x[z] (resp. y[z]) is the 
z-transform of (x(n)) n6Z (resp. (y(n)) neZ ). 

C. Synthesis FB 

The polyphase transfer matrix of the synthesis FB: H[z] = J2eH(ty z ~ i > satisfies: 

x[z] =H[z]y[z], (3) 

where the polyphase vector of the ouput signal of the synthesis FB (x(n)) nG z is defined similarly to 
(x(n)) ne g. We deduce from © that: 

M—l oo 

VnGZ,Vi€{0,...,JV-l}, x(nN-i)=J2 ^ H id {n- £)yj{l), (4) 

j=o e=-oo 

where H(£) = yHi t j(£)j ^ N0< - M ' Expressing (0]) with impulse responses, we can write: for every 
i G {0,...,N- 1} and n G Z, 

M-l oo 

x(nJV-i)=E E hjiNin-l)-^^) (5) 

j'=0 £=-oo 

which, by identifying (0]) and (f5]), allows us to deduce that 

H(£) = \hi(N£ - i)} , ^gZ. (6) 

w I Jy J iO<i<N,0<j<M 

These expressions hold for any oversampled FIR FB. 
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III. Inversion 

A. Invertibility of an analysis FB 

This work being focused on the construction of FIR synthesis filters, a preliminary point is the 
verification of the given analysis FB FIR invertibility. The polyphase representation of FBs offers the 
advantage of relating the perfect reconstruction property to the invertibility of the polyphase transfer 
matrix Il42l . The latter matrix belongs to the ring C[z, z~ 1 ] MxN of Laurent polynomial matrices of 
dimensions M x N. We emphasize that we do not look for any inverse MIMO filter, but for an inverse 
polynomial matrix in C[z, z~ l ] NxAI instead. In other words, we aim at obtaining a (non-necessarily 
causal) FIR synthesis FB. 

A first answer to this FIR invertibility problem can be conveyed through the study of the Smith 
McMillan form of a polynomial matrix ll43l . ll42l . but unfortunately this decomposition is quite costly. 
Park, Kalker and Vetterli also devised a method using Grobner bases |28l to study the invertibility 
of polynomial matrices which is applicable to the general multidimensional case. We describe here an 
alternative cost-effective method in the one-dimensional case. The following result gives a necessary and 
sufficient condition for a matrix to be left invertible, and thus, for the existence of such an inverse system: 
let H[z] G C[z, z ~ 1 j MxN be a polynomial matrix with M > N. The following conditions are equivalent: 

1) H[z] is "coprime", which means that the determinants of the maximum minors (sub-matrices of 
size N x N) are mutually relatively prime. 

2) H[z] is left invertible in the sense that there exists H[z] 6 C[z, z ~ 1 ] NxM such that = /jy. 
A proof of this result can be found in ll35l for instance. 

The first condition is directly applicable in practice to resolve the left invertibility of the polyphase 
transfer matrix. Using the following procedure, we can check numerically whether this condition is 
satisfied: 

® Extract a maximal sub-matrix H e [z] of H[z\. 

® Compute det (H e [z]), and determine its set of roots S e . 

© Consider another maximal sub-matrix. Remove from S e the elements which are not roots of the 

determinant of this sub-matrix. 
© Repeat step © until S e = or all maximal sub-matrices have been extracted. 
© If S e = then the polyphase transfer matrix is left-invertible; otherwise, it is not. 
The corresponding algorithm is easily implemented, leading to extract the roots of a single polynomial 
and check the roots of at most (^) — 1 = N \^L N y — 1 polynomials. If the polyphase matrix is left 
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invertible, the number of considered polynomials in practice is usually much smaller than ( N ) — 1, this 
bound being reached only when the matrix is not invertible. Note that in the case of causal filters (i.e. 
both H[z] and H[z] are polynomial matrices in C[z~ 1 ] NxM ), simpler invertibility conditions exist by 
invoking the so-called column-reduced property fi4l . B31 . Also notice that, one of the advantages of 
this algorithm over other methods is that it can be fully numerically implemented. 



B. Computation of an inverse FB 

The method proposed in Section IIII-AI only guarantees the existence of a left-inverse, corresponding 
to an FIR synthesis FB. Since it does not provide a constructive expression, we now perform the actual 
computation of an inverse polyphase transfer matrix. We assume hereafter that H[z] was proven to be 
FIR left invertible. 

Since the goal is to achieve PR, we search for a matrix H[z] in C[z, z ~ l j NxM such that H[z) H[z] = In 
and there exists {pi,P2) £ N 2 such that the polyphase transfer function of the synthesis FB reads: 
H[z] = J2eL- Pl H(£)z~ e . The resulting overlapping factor of the synthesis filters is p = p± + P2 + 1. 
When working with Laurent polynomial matrices, these integers p\ and p2 are a priori unknown, whereas 
with polynomial matrices a bound exists ll44l . By rewriting the PR property in block convolutional form, 
we get the following linear system: 

HH = U 



where 



On, Pi n, In, 0jv",(p 2 +fc-i)iv 
/ 

and H T = 



H(0) 



g K iVx(fe+p-l)AT yJ = 

• H(k-l) 



H(-pi 
\ 



,H( P2 ) 



£ C 



NxpM 



£^pMx(k+p-l)N 



(V) 



(8) 



(9) 



y H(0) ••• H{k-1) J 

As already mentioned, p\ and p 2 are unknown, but since the system $7) is supposed invertible, at least 
a couple of integers (pi,P2) solving the system exists. The values of p\ and p2 are actually obtained 
by increasing the value of p and looking for every couple satisfying p = p\ + p2 + 1, starting with 
p = 1. Hence, for a given p, we consider all (pi,P2) in {(p — 1,0), (p — 2, 1), (0,p — 1)}. The first 
p allowing a Moore-Penrose pseudo-inverse ||46l solution to © provides an inverse polyphase transfer 
matrix of minimum order. 
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C. Hermitian symmetric case 

1) Symmetry conditions: It is well known that the Fourier transform of a real signal is Hermitian 
Symmetric (HS): its frequency decomposition is symmetric for the real part and anti-symmetric for the 
imaginary part. Conversely, if the coefficients are HS in the frequency domain, then the reconstructed 
signal is real. This property is very useful for real data filtering, which often consists of removing or 
thresholding coefficients in the frequency domain before reconstructing. Securing the reconstruction of 
real-valued signals from the transformed coefficients is thus a desirable property. In this section, we study 
the HS case and its effects on the methods proposed in the previous sections. 

The HS property in the synthesis filters is satisfied provided that, considering any symmetric subband 
indices jt G {0, M — 1} and M — 1 — jf, for any coefficients (yi(n))o<i<Af such that yi(n) = if 



(i,n) / (jf,n f ) or (i,n) / (M - 1 - jf,n f ) with n f G Z and, such that yj f {n f ) = y M -i-j f {nf), a 
real-valued signal is reconstructed. The reconstructed signal reads: 

M-l oo 

x(m) = Y M m ~ N t)Vj(£) = h jf (m - n f N)y jf (nf) + h M -i-j f (m - n f N)y jf (n f ). 

j=0 £=-oo 



A necessary and sufficient condition for x(m) G R for all yj f {nf) G C, is that hj f ( 



m — n 



h 



M-l-j 



Am — rifN). This condition must be verified for any couple of integers (jf,rif). The condition 



on the synthesis filter is then: Vj G {0, M — 1} and Vn G Z, hj(n) = fiM-i-j{n). Using ©, we 
rewrite the condition as 



G {-pi, ...,p 2 } 



where J m is the M x M counter-identity matrix: 

/ 



' M 







1 



H(£) =H(£)J M - 



(10) 



1 







We have supposed here that the transformed coefficients of a real signal exhibit the HS property. In other 
words, for any real signal (a;(n)) ne z, and for any couple (jf,rif) G {0, ...,M — 1} x Z, the output of 



the analysis FB verifies: yj f (rif) = yM-j f -i(nf). This condition can be rewritten: 

hj f (m)x(Nnf — m) = hM-j f -i(,Tn)x(Nnf — m). 



Considering a zero input signal x except for one sample, we deduce that hj. (n) = h 



lM-jf-1 



(n), which 



is equivalent to 



H(£) =J M H{£). 



(11) 
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Hence, if the analysis FB verifies Condition (fTTT) . then the coefficients after decomposition satisfy the 
HS property. 

Remark 1: Consider an invertible HS analysis FB. By inserting (fTTb in the PR condition we get: 

min(p2,^) 

5tI N = £ H(s)J M H(£-s). 

s=m&x(£—k+l,—p 1 ) 

It implies that if H = H(—pi),--- ,H(p 2 ) is a solution of the linear system (0 then, under the 



is also a solution of the linear 



HS hypothesis on the analysis FB, TL2 = H{—Pi)Jmi ' ' ' ,H(p 2 )Jm 
system. Finally, it follows that the sum: Ho = \{H + H2) is also a solution. Moreover, this solution 
verifies Condition (fTOb by construction. 

In other words, we have proved that an invertible HS analysis FB admits at least one HS synthesis FB. 

2) Construction method: We suppose here that the analysis FB was proven invertible and that the 
matrices H{£) satisfy Condition (fTTT) . Our objective is to build a synthesis FB possessing both the PR 
and HS properties. 

a) First case: M is even: First, we rewrite Conditions ( fTOl) and (fTTT) : 



P^We {-pi,...,p 2 } 



H{£) = [H 1 (£),H 2 (£) 
H 1 (£)=H = 2(£)Jm' 



Hi e C 



NxM' 



and H2 E C 



NxM' 



TT]) ^ ye e {o,...,fc- 1},< 



H(£) 



Hi(£) 
H 2 (£) 



H! E 



and H 2 E 



^M'xN 



_ Hi(t) =J M ,H 2 (£) 
with M' = M/2. Combining these conditions and the PR property, we get: 

min(p 2 ,^) min(p 2 ,^) 

5 £ I N = E Hi(8)Hx(l-s)+ ]T H^Jm.H^-s). 

s=max(£— fc+l,— pi) s=max(^— fe+1,— pi) 

Since 7^, = 7ms the previous equation can be seen as the sum of a complex matrix with its conjugate, 
leading to a real matrix. We deduce that 



(12) 



'-8,1 



min(p2,^) 

E 

=max(f — fc+1,— pi) 



H?(s),-H{(s) 



( Hf(£-s) 



where A R is the matrix of the real part of a matrix A and A 7 is its imaginary part. We will then define 
the following matrices: 



H 



H?(- Pl ), -H[(- Pl ), ■ ■ ■ ,tff (J*), -H[{ P 2) 



iNxpM 



(13) 
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and 



( fff(O) ••• H?(k-l) 
H{(0) ••• tf{(fc-l) 








fff (0) • • • Hf{k - 1) 

H{(0) ••• H{(k-1) ) 

b) Second case: M is odd : Similarly to the first case, Conditions (fTOl ) and (fTTT ) can be rewritten: 



V 



x(fc+pi+p 2 )7V 



ToJ^wg {-pi,...,p 2 }, 



H( 



H 1 (£),c 1 (£),H 2 (£)\ , H 1 G C JVxJkf ' and tf 2 G C JVxM ', 



= H 2 (£)J M , and ci(^) G 



flU OWG {0,...,fc- 1}, 



C2^ T 



#1 G 



and #2 G 



iM'xN 



, =Jm>H 2 {£) and c 2 (£) G 1^ 

with M' = (M — l)/2. Combining these conditions with the PR equation and following the same 
reasoning as in the previous section, we deduce that 



1 



min(p 2 ,£) 

:5pIn = ^2 

s=max(£— fc+1,— pi) 



^ H?(£ - s) ^ 

c 2 (l-s) T 



Subsequently, we introduce in this case the following matrices: 



flf(-i>o. 



jJVxpM 



(14) 



and 



/ fff (0) 

C2(0) T 

fff(0) 



(Jfe - 1) 

c 2 (fc-l) T 

H[{k-1) 






V 



tff(o) 

C2(0) T 

H{(0) 







Hf{k - 1) 

C 2 (fc-1) T 

H[(k-l) ) 



x(k+p 1 +p 2 )N 
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c) Conclusion : In both even and odd options, we solve a linear system of the same size as the one 
of Section UlI-BI but with real coefficients in this case. More precisely, with the introduced notations, we 
have: 



TL S TL S =U S = - 



0n, Pi n, In, N ^ P2+k _i) N 



The system is then solved, in the same way as in Section IIII-BI For increasing values of p (starting with 
p = 1), for each couple (pi,P2) G N 2 such that p = p\ + P2 + 1 we try to invert the generated system 
through a Moore-Penrose pseudo-inversion. 

IV. Optimization 

A. Dimension reduction 

1 ) General case: Before addressing the issue of optimization in itself, let us rewrite the linear system 
expressing the PR property. The analysis FB is still supposed invertible. Let r be the rank of the matrix 
TL £ £(k+ Pl +p 2 )Nx P M _ Wfi assume( j t hat r < Mp (with p = p x + p 2 + 1). Performing a Singular Value 
Decomposition |47l (SVD) on this matrix yields TL = UqI^qVq , where So G C rxr is an invertible 
diagonal matrix, U G C Ar(7c+p_1 ) xr and Vb G C Mpxr are semi-unitary matrices (i.e. UqU = I r and 
y*y = I r )- Therefore, there exists U\ G C N(k+ P -i)x(N(k+ P -i)-r) and Vi £ C M P x(M P -r) such that 

\Uq, U\\ and [Vo, Vi] are unitary matrices. When an inverse polyphase transfer matrix exists, a particular 
solution to (0 is 7i° = H}lA , where li} = VoSq Uq is the pseudo-inverse matrix of TC. Equation $7$ 
is then equivalent to Uq'EqVq {TL — H ) = 07v(fc+p-i)xAf- Since UqUq = I r and So is invertible, we get: 
Vq (H - H°) = rxN . In other words, the columns of H-H° belong to Ker(Vo), the null space of Vq. 
Moreover, it can be easily seen that Ker(Vo) is equal to Im(Vi). We then obtain the following affine 
form for TL: 

H = V 1 C + fL°, (15) 

where C G C ( - m p~^ xN . 

The construction of a synthesis FB thus amounts to the choice of C. If C = 0(M P -r)xN> then the 
obtained synthesis FB is the PI FB. This expression can be further rewritten into a more convenient form 
for optimization purposes. First, we define the matrices (Vj)j g /o j ... ) M-i} by: for all I G {—pi, ■ ■ ■ ,p 2 } and 
n G {0, . . . ,Mp-r-l}, (Vj)e +Pun = Vi((£+pi)M +j,n) , with Vi = [Vi(s, n)] < s< Mp,o<n<M P -r- Ac- 
cording to ([151) and dHJ), we can write for all £ G {—pi, ■■■,P2}, i G {0, ...,N - 1} and j G {0, M - 1}: 

Mp-r— 1 

HiM= E (^W,nC(M+^.(£), 
n=0 
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where (Hij(£)) _ p <t<p represent the impulse responses of the synthesis FB, (H^iJPj) <i< p corre " 
spond to the PI solution and C = [C(n, i)]o<n<Mp-r,o<i<N ■ For all j G {0, M — 1}, we introduce the 

Hfjil) for all I G {~Pi,...,p 2 } and i G {0,...,N- 1}. We thus 



matrices Hj defined by: \Hj 
obtain 



£+Pi,i 



H i , j {i) = (V 3 C+H 



(16) 

This equation is used in Section IIV-B 1 1 to simplify the optimization problem raised by the design of the 
synthesis FB. 

The above expressions are given in the complex case, but they naturally remain valid in the real case. 
This will be illustrated by the first example of Section IV-B3al 

2) Symmetric case: In this section, we adapt the results of the previous section to the HS FB case. 
The notations used here are similar to those introduced in Section IIII-C2I It is worth noticing that we 
can calculate the matrix Tt directly from 7i s , as defined in Section ITII-C2I when M is either even or odd, 
in the following way: 

H = P rc H s (17) 

where the matrix P rc G C pMxpM is the block diagonal matrix built with the block: 

( I I 

1M> —Wm> 
Jm< lJ M> 

if M = 2M' (even case, as seen in Section IIII-C2ab and: 



V 



( I M > -%I U , ^ 



V2 







V 





Jw iJm 1 J 

if M = 2M' + 1 (odd case, as seen in Section [Tll-C2bb . By applying once again an SVD on TL S , and by 
following the same steps as in Section HV-AU we end up with an equation similar to ( fT5l ): 



H s = ViC + H s 



(18) 



Note that, according to the properties of the SVD, the matrix C is now real-valued. By noticing that 
Tt° = Pvc^ri^ and setting Wi = PrcVi, we finally obtain: 

n = W 1 C + n°. (19) 

We next define the matrices (Wj)o<j<M-i- for all I G {—pi, ■ ■ ■ ,p 2 } and n G {0, . . . , Mp — r — 1}: 
(Wj)i+p 1>n = Wi({£ + pi)M + j,n). Using £TJ]> as in Section IIV-A1I we get: 



ffy-W = [WjC+H 



£+Pi,i 



(20) 
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B. Optimal solution 

1) General form for the cost functions: Depending on the desired properties for the synthesis FB, 
several cost functions can be employed. We first propose a generic cost function formulation and then 
provide practical examples based on the filter time or frequency spread, respectively. 

Our goal is to optimize the filter shape given by the coefficients h of the synthesis FB, subject to 
the perfect reconstruction property. According to the results in Section IIV-AI it is possible to represent 
the coefficients in the general case by using ( fl6l ). The optimization favorably takes place in the reduced 
dimension space the matrix C belongs to (compared with the dimension of the space of the coefficients 
of h), thus allowing us to reformulate the optimization under a perfect reconstruction constraint as an 
unconstrained problem. In this context, the generic cost function form we consider is: 



J(h) = J(C) = £ 



M-i \\vfi+H] 



2 



j=o WjC+H] 



Hereabove, the following notation has been employed: 



VAeC" x * \\A\\ 2 K = J2 AijA^eK&i', £,£'), 

(i,i',£,£') 

where K and A are (N xJV xpxp) kernels. Moreover, we assume here that \\A\\ K represents a semi-norm 
over C Nx p and it is thus real non-negative. Let K' be the matrice defined by ^-|_^jv,i'+^'iv = 

K{i,i',e,e) 

for all (£, £') G {0, ...,p — l} 2 and (i, i') £ {0, N — l} 2 . Without loss of generality, this matrix can be 



taken positive semi-definite, which implies that K' i+eNi , +e , N = K' i , +l , Ni+lN and, thus, K(i,i',£,£') 



K(i', i, £', £). We deduce the following expression: 

\\A\\\= A he A~^K(i,i', £,£')= A itt Ai lt t,K(i',i,£',i). 

,£>£') 

This relation will be used to simplify some equations in Section IIV-CI and Appendix We finally notice 
that K' is a positive definite Hermitian matrix if and only if is a norm. 

2) Impulse responses optimization: A first objective is to obtain impulse responses (^j)o<j<m f° r tne 
synthesis filters well-localized, around some time-indices (rfij)o<j<M- We now explain the link between 
the cost function form introduced in the previous section and the previously described dimension reduction 
to further simplify the problem. 

The considered cost function is the following: 

Em \ m ~ rnj\ a hj(m) 

M h > = 2^ Wt J ~ 

i=o J2 m hj(m) 
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with a£i; and weights (wy)o<j<j\/ G (R + ) M such that Ejii^tj = 1. 



If a = 2 and 



E m m 


| ( m " 




Ej 


^3 ( m ) | 





, then J\(h) represents a weighted sum of the standard temporal 



dispersions measuring the time localization of a filter hj ll48l . Combined with Q, we get: 



M h ) = E ^tj 

3=0 



^ E?=_ Pl Eilo 1 l^r - i - m 3 r H id (£) 
^=0 E£_ Pl ES 1 

We now introduce the kernels K { - and A defined by 

K}(i, i',£ + Pl , £' + pi) = |&V - i - r%r 
A(i,i',£ +pi,£' +pi) = Si-i'St-t', 
for all j S {0, ... ,M- 1}, E . . . ,p 2 } 2 and (»,*') £ {0, . . . , iV - l} 2 . Using ([111) we write: 



(21) 



P2 JV-1 

E E 

!=-pi i=0 



and 



Hi A*) 



P2 N-l 

^ E E( £Ar -* 

Hereabove, ||.|| A reduces to the Frobenius norm. Finally, we deduce that 



M-i \V j C+H j 
Mh) = E 



2 

K 1 



i=o VX+H, 



MC). 



The constrained minimization of J\ is then reexpressed as the unconstrained minimization of J\. 

3 ) Frequency response optimization: We proceed similarly to the previous section, for a different cost 
function J\(h). Our goal, dual to that in the previous section, is now to regularize the frequency responses 
of the synthesis FB by concentrating the frequency response of each filter hj around some frequency fj. 
This is achieved by minimizing 



Af-l J 
Mh) = E 

j=0 



-1/2+/, \ v h 



h 3 \v] 



dv 



'1/2+ fi 
1/2+A 



(22) 



where a £ R^_, (a»fj)o<j<Af G (K + ) M with Y^fLo 1 ^t,j = 1 an d, hj[.] is the frequency response of the 
j m synthesis filter, defined as 

P2 N-l 

v»e [-i/2,i/2[, *,-M = E E ^i,iW e_2,7r(w " i)l '' 

£=-pi j=0 
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When f j 



-1/2 
-1/2 



hj[v] 



dv 



1/2 
/2 



hj[u] 



dv 



-, the cost function J\{h) represents a classical weighted frequency disper- 



sion measure for the synthesis filters. We then define the kernel 

"1/2+/, 



' J -1/2+ f \ 



Hi 



1/2+/, 
1/2 



-1/2 

i2 w r t2 



I/ |« e -2 M r(JV(<-/')-(i-i'))(«'+/ i ()j I/j 



(23) 



with (i,i',£,£') e {0,...,iV- l} 2 x {- Pl ,..., P2 }' 

Remark 2: The examples provided in Section IV-B3I are obtained with a = 2. In this case, the explicit 
expression of the kernel becomes: 



12 

Wfi (_l)A r (^')-(*-i') e -2* 7 r(iV(£-£')-(i-i'))/ J 



2tt 2 (N(1-1>) - (t-*')) 2 

with e {o,...,iV- i} 2 x {-pi,..., P2 } 2 . 

Combining these notations and ( fT6l ), we have 

2 



if i = i' and 
otherwise 



'•1/2+/, 
-1/2+/, 



V<C+H, 



Kl 



Invoking Plancherel's theorem and the kernel A defined in Section IIV-B21 we obtain: 

P2 N-1 



/•1/24 


-/» 




2 










1-1/2- 









«=-p! i=0 



Finally, substituting these expressions in (1221 ) yields 

W = E 



A'' 



Once again, the constrained optimization problem has been reformulated as an unconstrained one. 



C. Gradient optimization 

The constrained optimization problem being turned into an unconstrained minimization, we now provide 
more details about the minimization algorithm we employ. In this work we have used a simple gradient 
algorithm with an adaptive step \x n . The algorithm can be summarized as follows: 
® Initialization: Co = 0, n = 0. 

® Mn = 1 

® Computation of D n = VJ(C re ). 
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© While J(C n - p n D n ) > J(C n ), set p n <- -X^. 

© C n +i = C n — [X-nPn- 

® If ||C n +i — C n \\ > e then increment n and go to step ®. 

The step-size p n used here remains large as long as the algorithm is getting closer to a local minimum 
(in other words, as long as J(C n+ \) < J(C n )). It is only adapted (reduced) to prevent the criterion from 
increasing. The initialization with Co = entails that we consider the pseudo-inverse synthesis FB as the 
starting point for the algorithm. In practice, e was set to 10 -13 . 

Other step selection strategies exist: constant or optimal steps, steps satisfying Wolfe or Armijo 
conditions |49l , ll50l . The method used in this work is easy to implement and is well-suited to the 
different cost functions we have considered, while keeping a reasonable complexity. 

As the cost functions considered in this work are not convex, there is no theoretical guarantee that 
the algorithm converges to a global minimum. Yet, as is shown in Section [V] initializing this method 
with the PI synthesis FB provides quite good results and extensive simulations have confirmed this good 
behaviour. 

The expression of the gradient for the general cost function is given in Appendix |A] and is next applied 
to J t an d J\ in Appendix |B] 



D. Optimal solution in the symmetric case 

1) Cost functions: Using the same notations as in Section UV-B1 in the HS case, the following form 
of the cost function is found: 



M-i WWX+H, 



j s (h) = j s (c) = J2 

7=0 pTjC+ffj 

As in the general case, (l20b has been used to transform the constrained optimization problem on h into 
an unconstrained minimization problem on C. 

2) Examples of cost functions: Equation (l20b is very similar to (fT6l ). We consequently define the 
cost functions in the HS case following the same approach as in Sections IIV-B2I and IIV-B3I Thus, the 
following functions are considered: 

i0 



m-i WX+H, 



3=0 WjC+ff, 
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to concentrate the time localization of impulse responses, and: 



M i \\Wf + H j 



.0 2 



i k< 

3 



3=0 WjC+Hj 




to enhance the frequency selectivity of the filters. Their gradients are provided in Appendix O 

V. Examples 

As emphasized in the introduction, a wide variety of filter banks and design choices can be made. In 
this section, we have chosen to work with three different examples exhibiting interesting properties and 
allowing us to show the benefits incurred in the proposed inversion and optimization methods. 

A. Considered filter banks 

1) Real Lapped Transforms: The study, developed for the general complex case, remains fully ap- 
plicable to the design of real filter banks. As an illustration, we first consider real lapped transforms 
introduced in the middle of the 90s under the name of GenLOT (generalized linear-phase lapped orthog- 
onal transform) |5TTl . Those transforms generalize the DCT (Discrete Cosine Transform) and the LOT 
(lapped orthogonal transform). 

To illustrate the inversion method, we have chosen a GenLOT with M = 16 filters of 32 coefficients. 
This FB is invertible, in a non standard oversampled use, with parameters N = 8, k = 4 and k' = 2. Its 
impulse and frequency responses are represented in Figure [2 This FB is real and does not satisfy the HS 
condition. By using the method described in Section IIII-BL we find p\ = 3 and p2 = (hence p = 4). 
The frequency and impulse responses of the synthesis FB computed with the pseudo-inverse are shown 
on Figures Oa) and @Ia). 

2 ) Modulated complex lapped transform: We now consider another analysis FB based on a windowed 
generalized Fourier transform, corresponding to a modulated complex lapped transform (MCLT). This 
family of FB has been used by Kingsbury 01 or Malvar Q for applications in video as well as audio 
processing. The analysis impulse responses are: hi(n) = E(i,n)h a (n), where 



and (/i a ( n ))i<n<fc7V is an analysis window. In this paper, we consider two analysis windows. The first, 



E(i,n) 



1 



^-^+|)(n-^+f)^ 



e 



defined by 



Vn G {1, kN} , h ai (n) = sin 



( 



kN + 1 



) 
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Fig. 2. (a) Impulse and (b) frequency responses of a GenLOT analysis FB. 



is a standard sine window, employed for example in |7J, J3J. The second (h a2 (n))i< n <kN, corresponds to 
a zero-phase low -pass filter with cutoff frequency 2ir/(kN), built from a Kaiser window. This window, 
with better tapering than h ai , was used for instance in |[32l . It is interesting to note that this analysis FB 
family, with both analysis windows, satisfies Condition (fTTT) . In other words, it can be used to illustrate 
our approach in the HS case. The method from Section IIII-AI was employed to verify the invertibility on 
this FB, with both analysis windows and parameters N = 8, k = 3 and k' = 7/4. We then compute a 
first synthesis FB with the PI method of Section IIII-BI For the analysis FB with h a± window, the minimal 
parameters pi = 2 and P2 = were obtained. The frequency response of this synthesis FB is represented 
in Figure [3a). In the h a2 case, the minimal parameters pi = 2 and P2 = were, once again, found when 
applying the method of Section IIII-Bi An HS synthesis FB is then derived from this filter bank using 



July 21, 2009 



DRAFT 



18 



4r 

2 ' 
- 

0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 



0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 



(a) 




0.1 0.2 0.3 0.4 



4r 

2 

- 

0.1 0.2 0.3 0.4 

2 
1 

2 
1 

2 
1 

2 
1 




0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 



0.1 0.2 0.3 0.4 



2r 
1 

r 

2 
1 


0.1 0.2 0.3 0.4 






0.1 0.2 0.3 0.4 



0.1 0.2 0.3 0.4 




0.1 0.2 0.3 0.4 



0.1 0.2 0.3 0.4 



Fig. 3. First example (with the GenLOT FB of Section |V-AT 
the pseudo-inverse method and (b) after optimization with cost 



(b) 

\: frequency response of the synthesis FB obtained (a) through 
function Jt. 



the method of Section ITII-C2I to directly build an HS synthesis FB. The frequency and impulse responses 
of the resulting synthesis FBs, in the h &2 case, are shown in Figures [6] and |7J respectively. Figure 0a) 
shows that synthesis filters present a symmetric behavior for their coefficients (in other words, they have 
a linear phase) while the synthesis FB in itself is not HS. We also notice that the frequency selectivity or 
time-frequency localization of the filters obtained through the pseudo-inverse methods is not satisfactory. 

B. Optimization examples 
1 ) Kernel parameters: 

a) Temporal kernel Kj: It is defined by (|2"T1 ). The parameters rnj define the (temporal) positions 
around which the impulse responses of the j^ 1 filter should be concentrated. To obtain well tapered 
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Fig. 4. First example: impulse response of the synthesis FB obtained (a) through the pseudo-inverse method and (b) after 
optimization with cost function J\. 



filters, we need to concentrate the impulse responses around the middle of the filter support. Therefore, 
the same parameter was used for all filters. The support of the filters being {— p\N — N + 1, ...,p2N}, 
we have chosen rrlj as 

wn . . . M _ _ P2N + 1 - piN - N 
VU < J < M, rrij = m = . 

In our design example, a has been set to 2. 

b) Frequency kernel K^: It is defined by (123T ). The parameters fj represent the reduced frequencies 
around which we want to concentrate the frequency responses of the synthesis filters. More precisely, 
we chose fj such that it is centered inside the bandwidth of the analysis filter hj. The exponent a has 
been set to 2. 
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c) Weight parameters: In the proposed cost functions, the parameters cjy and control the relative 
importance of the different filters in the optimization process. For the following examples we have chosen 
equal weights: 

Vj G {0,...,M-1}, wy =Wfj = — . 

In other words, we aim at obtaining synthesis filters with similar behavior. 

2) Computation time: In Section IIV-A1 we have seen how to parameterize the system and thus 
how to reduce the dimension of the optimization problem. To evaluate the gain resulting from this 
parameterization, Matlab programs were written to compare the solutions of the constrained problem, 
using function fmincon, with the solutions of the unconstrained problem, using function fminunc 
and using the gradient method explained in Section HV-CI The two functions fmincon and fminunc 
were chosen as examples of optimization implementation, while the gradient procedure can be easily 
applied in different languages without requiring Matlab. These programs were tested with the analysis 
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Fig. 6. Modulus of the impulse responses of the synthesis FBs (in the MCLT case with h a2 window): (a) pseudo-inverse and 
(b) symmetric version with the method of Section HII-C2I 



FB introduced in Section IV-A2I with the h ai window and the following parameters: overlap factor k = 3, 
redundancy k' = 7/4 and downsampling N £ {4,8,16}. The cost function used was J t (as defined 
in Section IIV-B2I ). Table U shows the computation time for the different methods on a computer with 
2.16GHz Intel Core2 T7400 CPU and 2Gb of RAM. 

TABLE I 

Computation time to optimize a synthesis FB with different methods using Matlab. 





N = 4 


N = 8 


N = 16 


Constrained optimization (with fmincon) 
Unconstrained optimization (with fminunc) 
Unconstrained optimization (gradient algorithm) 


1.2s 
0.04s 
0.06s 


120s 
0.7s 
0.6s 


8800s 
8s 
7s 



A first interesting result is that all three methods, starting from the same FB, converge to almost 
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identical synthesis FBs. The computation times are however very different: more than two hours (with 
N = 16) for the constrained optimization against a few seconds for the unconstrained optimizations. We 
can also notice that the gradient algorithm is as fast as the fminunc Matlab function. This shows that 
the optimization method can be easily implemented, through a gradient algorithm, with no performance 
loss and without having to resort to the fminunc Matlab function^. In other words, this last result 
indicates that in some applicative contexts in which Matlab is not available the optimization method can 
still be easily and efficiently implemented. 

3 ) Examples of optimized FBs: In this section, we present optimization result obtained with the 
different FBs introduced in Section |V-A| and using the different proposed cost functions. 

'This fairly sophisticated function uses an interior-reflective Newton method 1521 . 

2 A matlab toolbox for FB optimization is available here: http://www.laurent-duval.eu/misc-research-codes.html 



July 21, 2009 



DRAFT 



23 



a) General case: We have applied the optimization method on the real FB introduced in Sec- 
tion [VjAl] with parameters: N = 8, k = 4 and k' = 2. The employed cost function is J\. The result is 
shown in Figure Hfb), the coefficients before optimization (obtained with the pseudo-inverse method) are 
also displayed. It is clear that the impulse responses of optimized filters are better concentrated around 
the middle of the support. Figure [3] illustrates that the gain in time localization does not entail a too 
severe loss in frequency selectivity of the optimized filters. 

A second optimization example is given using an MCLT FB with analysis window h ai and parameters 
N = 8, k = 3 and k' = 7/4. The resulting frequency responses after optimization with cost function Jf 
are represented in Figure [3b). We observe that the frequency responses after optimization exhibit more 
regularity and an improved selectivity. The proposed cost functions take into account all synthesis filters 
at once. It therefore interesting to look more closely at each filter independently and determine whether 
the optimization leads to better results. In Table [III the frequency dispersion of each filter is reported 
before and after optimization with cost function Jf. In this case, the overall frequency dispersion of the 
optimized filters has been noticeably improved and spread variability has been drastically reduced. 

TABLE II 

Frequency dispersion of the synthesis filters optimized with cost function Jf. 



Filter 


Freq. disp. 


Filter 


Freq. disp. 


h 
hi 
h 2 

h:> 
hi 
h 5 
h 6 
h 7 
h 8 
h<j 
hio 
hn 
hi2 
his 


0.0290 
0.0851 
0.0569 
0.0606 
0.0658 
0.0596 
0.0363 
0.0156 
0.0411 
0.0231 
0.0508 
0.0499 
0.0574 
0.0521 


h°r 
hi vt 

K pt 
K pt 
h°t 

ir 

h° vt 

h r 

,opt 
ll 10 

"ll 

Topi 
n 12 

Topt 
"13 


0.0111 
0.0110 
0.0109 
0.0112 
0.0110 
0.0109 
0.0111 
0.0111 
0.0109 
0.0110 
0.0112 
0.0110 
0.0109 
0.0111 


Sum 


0.6833 




0.1544 



July 21, 2009 



DRAFT 



24 



b) Results in the symmetric case: The optimization procedure was next applied in the HS case to 
the FB of Section IY-A2I (with analysis window h a2 ). In this case, the cost functions J\ s and J\ s were 
employed. Once again, the following parameters were used: N = 8, k = 3 and k f = 7/4. Figures [8] and 
[9] show the optimization results. We observe that the optimizations with these two cost functions lead 
to FBs with different characteristics: as expected, with J\ s the impulse responses are better concentrated 
than with J\ s and, conversely, with Jf S the frequency selectivity is better than with J ts . 
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Fig. 8. Third example (in the MCLT case with window h a2 ): modulus of the impulse responses of synthesis filters (a) before 
and after optimization (b) with the cost function Jt s and (c) with J fs . 



4) Comparison: To conclude this example section, we propose a comparison with an existing filter 
bank design. We have chosen to compare our design methodology with the filter banks used in (71, (H. 
In these works, the considered FBs correspond to a modulated complex lapped transform with overlap 
factor k = 2 and redundancy k' = 2. For this application, our choice of N = 8 results in the filters 
shown in ||4). The synthesis filter bank is then built with a method equivalent to the weighted overlap-add 
technique. We have applied the methods proposed in this work to compute an optimized synthesis filter 
bank using the cost function J ts . In Table JIIJ the time dispersion of each synthesis filter computed as 
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explained in (H (hj^j ^ lg and after optimization {hfj 15 are re P orte d as wei l as the value 
of the cost function J ts . The time dispersion was clearly reduced with the proposed method. 

VI. Conclusion 

In this paper, we have proposed a method to test that a given oversampled FIR analysis FB is FIR 
invertible and a method to compute an optimized inverse FB. The optimization was performed for a class 
of cost functions allowing either to emphasize the time localization or the frequency selectivity of the 
filters. By rewriting the system defining the synthesis FB, we were able to parameterize the synthesis 
filters for a given filter length. This parameterization was then used to convert the constrained optimal 
synthesis problem into an unconstrained one, which can be solved with a simple gradient algorithm. 
The FB considered here are one-dimensional; it would be interesting to study how the proposed methods 
could be extended to the multidimensional case. Another perspective could be to study the case of FBs 
admitting an IIR left inverse that can be approximated using a FIR FB with very long support. 
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TABLE III 

Time dispersion of the synthesis filters optimized with cost function J t 



Filter 


Time disp. 


Filter 


Time disp. 


h , his 
hi, hi4 
h 2 , hi3 
h 3 , hi2 

/l4, /in 

h$, hio 
he, hg 
h-t, ha 


21.25 
21.25 
21.25 
21.25 
21.25 
21.25 
21.25 
21.25 


, opt 
n ' 

K pt , 
K v \ 
~h°*\ 
K pt , 
h° pt , 
K v \ 

h° v \ 


,opt 
"■15 
Topt 
"-14 

Kg 

,opt 

"■u 


2.4152 
2.1064 
2.3424 
2.5516 
2.5516 
2.3424 
2.1064 
2.4152 


X 


340 




37.6626 



Appendix A 
Expression of the gradient 



In this first appendix, we study the gradient of J (as defined in Section IIV-B lb with respect to C. We 



first need to calculate the gradient of f(C) 

df df df 
+ 1-. 



V j C+H j 



The matrix C being complex, we have 



dC m ,n dC^ n dCm : n 

E {{Vj)i,mh- n {VjC + H°jh,e + Wj)^j£'-n(VjC + H*)^) Kj(i, i 

„•/ pt V / 



,n 



iA't 



From this result, we deduce that 

VJ(C)) = 2V ^'A^'^ V i C + *fy M K 3& ^ n ) ~ <*jA(i,i',£,n)) 

J m,n 

3=0 



(24) 



with a 



j 



and (3j 



VjC+H 
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Appendix B 
Examples of gradients 
A first example is the calculation of the gradient of the cost function J t . Applying (l24l ) to the kernels 

(Kj)o<j<M, we get: 

d f o 



r)C 

<-"-'m,,n 

= 2u %j (v*(T j Q(v j c+H° j ))) , 

where (Tj) it g +Pl = \£N — i — ffij\ a for all < i < N and — p\ < £ < P2- The symbol represents the 
Hadamard product (or pointwise matrix product). Finally we obtain: 



9 ^ YjiM^j ~ g>g Q (VjC + g,-)) 

VJt^J - ^L & : 

3=0 n 



where = 1 for all < i < N and < £ < p. 



The same study can be carried out for the cost function Jf, with kernels (Kj)o<j<M- Rewriting the 



result under a matrix form does not simplify the final expression in this case. Hence, the gradient reads: 



fvJ,(C)) -IV £.V.>(Vj)<-. m (VjC +Hj)i, f (0,A';(i,i',<,n) - a ,5^,6,^) 

V / m.n ' 

j=0 



0] 



Appendix C 
Gradient functions in the HS case 

Similarly to Appendix lAl we first compute the gradient of f(C) 
matrix C, the only difference being that the matrix C is now real: 
df_ = 



WjC + Hj 



with respect to the 



8Cr, 



By using this expression and the relation Kj(i,i' ,£,£') = Kj(i' ,i,£' ,£), we deduce the gradient of the 
cost function J s : 



(V Js(C)) = 2Re r£ El ' l '' £ ' {W ^'- JW -' C ± H ^'' e ' (tWh ^ gj ~ a i K ^ w ' gjj 



(25) 



with a 



j 



WhC+H, 



and (3j = WjC+Hj , for all < j < M. Using (|25T) . the calculation of 



the gradient of J\ s yields 



VUC) = 2Re I £ m^£j ~ ® + »7» 



3=0 
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where Tj is defined as in Appendix |Bj For the second cost function we find 




) 
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